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Abstract 

The steady-state simplified Pn (SPn) approximations to the linear Boltzmann equa- 
tion have been proven to be asymptotically higher-order corrections to the diffusion equa- 
tion in certain physical systems. In this paper, we present an asymptotic analysis for the 
time- dependent simplified Pn equations up to N = 3. Additionally, SPn equations of 
arbitrary order are derived in an ad hoc way. The resulting SPn equations are hyperbolic 
and differ from those investigated in a previous work by some of the authors. In two 
space dimensions, numerical calculations for the Pn and SPn equations are performed. 
We simulate neutron distributions of a moving rod and present results for a benchmark 
problem, known as the checkerboard problem. The SPn equations are demonstrated 
to yield significantly more accurate results than diffusion approximations. In addition, 
for sufficiently low values of N, they are shown to be more efficient than Pn models of 
comparable cost. 

1 Introduction 

The mathematical equation describing linear particle transport problems is the linear Boltz- 
mann equation. Its large dimensionality and integro-differential structure make this equation 
difficult to solve analytically and numerically. Since analytic solutions of the Boltzmann 
equation can be constructed only for simple geometries, a significant amount of effort has 
been invested to develop approximations and hence, calculate numerical solutions for real- 
istic, multidimensional problems. The spherical harmonics (Pn) equations are a standard 
approximation, already known from the beginning of the 20th century (first ideas for ana- 
lytical solutions in and attempts with focus on numerical computations in Ej)- A 
drawback of the Pn equations in 3-D, which has made them unattractive for practical appli- 
cations, is their complicated coupling and the large number of equations, growing as (N + 1) 2 . 
In the view of the low computational resources in the 1950s, it was inevitable to come up 
with simpler equations for the solution of realistic problems: Gelbard |15} IT6t [T7] therefore 
proposed the steady-state simplified Pn (SPn) equations, which were easier to implement 



(because they could be written as a system of diffusion equations) and whose system size 
increases in general geometries only linearly as (N + 1) (versus quadratically as the Pjy equa- 
tions do). However, Gelbard's derivation in 3-D geometries was purely ad hoc (by taking 
the 1-D Pjsf equations and replacing the 1-D spatial operators by its 3-D generalizations, i.e., 
gradients and divergence operators). Due to the lack of a theoretical foundation, the SPn 
equations were not accepted as an approximation to the transport equation until first theo- 
retical justifications were presented (asymptotic and variational analysis in [33J I2H1 EH E])- 
In the framework of a Galerkin finite element method, the well-posedness of the steady-state 
SPn equations is shown in [JT] for N = 1,3, 5, 7, where a proof of existence and uniqueness 
is provided. A detailed review of the SPn equations can be found in [31J . 

Originally intended for applications in nuclear engineering, the SPn equations are, in- 
deed, implemented and used for neutron transport problems nowadays [21 EJ [18, 9|. After 
first theoretical foundations for the SPn method were provided, a wide range of additional 
applications has been developed mainly in the past decade, e.g., radiative cooling of glass 
[301 [TT] , radiative transfer in tissue [22] , fluorescence tomography [23] , design of combustion 
chambers for gas turbines [351 136] . crystal growth of semitransparent materials [Tj, and photon 
and electron radiotherapy [24] . 

The majority of previous investigations has focused on steady-state SPn equations. One of 
the ideas for deriving SPn equations is to explicitly solve for odd moments and introduce them 
into equations with even moments. This is the reason why the steady-state SPn equations 
reduce to a hierarchy of diffusion equations. Keeping the time-derivatives in the odd moment 
equations, this procedure cannot be applied in the same way for the time- dependent case. 

To the authors' knowledge, the first formal asymptotic derivation for time- dependent SPn 
equations is developed in |10] , and Finite Element numerical solutions of this system are 
computed in [12] . An alternative strategy for the derivation of moment methods for the time- 
dependent radiative transfer equation is the method of optimal prediction [371 03] ■ It turns 
out that this formalism yields existing moment models such as Pn and diffusion closures. 
Additionally, it is shown in [37] that this approach can be used to derive variations of the 
parabolic SPn equations from [TP] , 

In this paper, we present an asymptotic analysis for time- dependent SPn equations, and 
we also explain how these equations can be derived in an ad hoc way. Our analysis makes use 
of a scaling which is different from [10] and leads to final equations which are not equivalent 
to those investigated in |10j . We want to highlight the differences and similarities between 
the approach therein and the work presented here: 

• Guided by the fact that, in steady-state, SPn approximations are diffusion equations, 
the authors in [10] derive time- dependent SPn equations which are parabolic PDEs. We 
follow a different approach in our asymptotic analysis and derive a system of hyperbolic 
PDEs for the time-dependent SPn equations. 

• The analysis in |10j is performed by a parabolic scaling in which the time-derivative 
is scaled by e 2 . As the final SPn equations are only accurate for e ~ 0, this asser- 
tion implies that temporal changes of the solution should be very small. In contrast, 
the asymptotic theory in this paper does not require a scaling of the time derivative 
operator. Hence, solutions of problems with large time-derivatives should also be ac- 
curate. However, both asymptotic theories assume that space-derivatives (scaled by e) 
are small. 
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• The derivation in [10J unfolds an ambiguity of how to define the 4>2 unknown. This 
ambiguity is only partly captured by the introduction of a free parameter a in the 
approximate system. Although leading to more flexibility, different choices of this pa- 
rameter a can only give results which differ in the magnitude of C(e 6 ). Moreover, in 
order to obtain a well-posed system a must be bounded by < a < 0.9. It turns out 
that approaching the lower or upper bound, numerical solutions of the regarding sys- 
tem diverge from the true solution and develop spurious shapes. Since it is not obvious 
which value of a is "optimal" in some sense this issue remains unclear for practitioners. 

The analysis discussed here is similar to the asymptotic derivation of the steady-state 
SPn equations. In this new analysis, there is less ambiguity and no free parameter. 

• The asymptotics in this paper and in [10J are performed only up to SP3. However, as 
it is additionally presented here how the ad hoc derivation by Gelbard |15|. [TBI [T7] can 
be generalized to yield exactly the same SPn equations as asymptotically derived up 
to N = 3, it is straightforward to obtain higher-order SPn approximations. 

The first asymptotic derivation of the time-dependent SPn equations was not developed 
until 2007 [10]. However, simulations with time-dependent SPn equations had been per- 
formed before. In practice, to keep the structure of diffusion equations for the time- dependent 
generalization, some simplifications were proposed or time-derivatives were dropped in cer- 
tain equations. A different possibility, which is already included in some codes nowadays 
(e.g., PARCS [9]), was to simply add the partial time derivative to each of the steady-state 
equations in first-order form. It will be shown in Section [2] that the time-dependent SP3 equa- 
tions, gained heuristically in this way, can be obtained from those derived in this paper by a 
similarity transformation. However, time-dependent SP3 equations are also used in practice 
(e.g., |18j ) which are not equivalent to the asymptotically derived equations from Section [3j 

This paper is organized as follows: In Section [2| a derivation of the time-dependent 3-D 
SP3 equations is presented which basically follows the lines of the ad hoc derivation given by 
Gelbard [151 EH Ej • I n particular, this approach is purely ad hoc. An asymptotic theoretical 
foundation is then given in Section [3] for the time-dependent method up to order N = 3. 
It is important to stress that our asymptotic analysis in Section [3] is only strictly valid for 
a homogeneous medium. In Section [2j heterogeneous media are considered in an ad hoc 
derivation. This procedure can be generalized to SPn equations of arbitrary order, and SPn 
solutions are compared numerically to diffusion and Pn results in Section |4j A discussion 
and conclusions are given in Section [5} 

2 Classic (Ad Hoc) Derivation 

In 1-D slab geometry, the Pn equations have a simple structure, and their number of un- 
knowns is only (N + 1). However, extending them to multi-dimensional geometries implies 
an expansion of the angular flux in spherical harmonics. Many extra degrees of freedom are 
added and an additional coupling occurs. Consequently, their original simplicity is lost and 
the number of equations increases quadratically in N. To preserve the pleasing form of the 
1-D slab geometry case, one can formally replace certain terms in the Pn equations in a 
certain way to obtain the simplified Pn approximation. 

The time-dependent, monoenergetic, isotropically scattering linear Boltzmann equation 
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reads as follows: 
Id* 



, y x,n,t) + n-v^(x,n,t) + E t (x,t)^(x,n,t) 

v at 



S s (x, t) 



4tt 



( ^{xM,t)d9! + ^Q(x,t), (1) 

J S 2 47T 



where the parameters are 



T, s (x,t) = scattering cross section, 
T, a (x,t) = absorption cross section, 

Ti t {x, t) = T, s (x,t) + S a (x, t) = total cross section, 

Q(x,t) = internal source. 

The angular flux *(x, Q,t) describes the particle density at position x E M 3 and time t 
traveling in direction G S 2 at velocity v. Penetrating the background medium, particles 
interact with the material, which is specified by scattering and absorption cross sections 
t) and T, a (x,t). However, the particles are assumed not to interact with each other. 
The first term of eq. ([I]) is the temporal rate of change in *, the second is the leakage or 
drift term, and the third quantifies the loss of particles due to out-scattering and absorption 
by the medium. The right-hand side of eq. characterizes the gain in particles: isotropic 
in-scattering processes are modeled by Y^ s (x,t) times the integral of \& over the unit sphere 
(all possible outgoing directions). Additional particles can also be inserted into the system 
by an internal isotropic source Q(x, t). 

To derive the SPn approximation to the above Boltzmann equation by a classic (ad hoc) 
procedure, eq. ([I]) is first restricted to planar geometry: 

1 <9* 

v at ox 

E s(x, t) f* [ ^( X:fi ' iip ' :t )d^dip' + ^-Q{x,t). (2) 



4vr Jq J_i 4ir 
Operating on eq. ^ by Jq W ■ dip and defining the azimuthally-integrated angular flux 

ip(x,fi,t):= ty(x,n,(p,t)d<p, 



JO 

we obtain the 1-D azimuthally-symmetric transport equation: 
1 dip dtp 

v~dj( X,fJ ' ,t ' +/i a^^ x '^'^ + ^t(x,t)ip(x,fi,t) 

= s^t) £^y^ )d//+ i Q(M) . (3) 

Multiplying eq. ^ by Legendre polynomials P n {^) and integrating this equation over — 1 < 
/j, < 1, we obtain, after using the recursion relation for P n Qu), 

1 d(j) n , . d ( n + 1 , , . n , . A . ,,, . 
~(x,t) + — - — —(/) n+ i(x,t) + - — — n _i(x,i) + T, t (x,t)<p n (x,t) 



v dt y ' ' dx V2n + 1 2n + 1 

= 8 n n (S s (x, i)0o(» } *) + Q(x, *)) , (4) 
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where we defined the Legendre moments of ip(x, fx, t) as: 

»i 



4>n(x,t) 



P n {ix)ip{x, fj,, t)dfj,, n > 0. 



(5) 



i 



Eq. Q is exact for all integers n > 0. Unfortunately, it never yields a closed system of 
equations; there is always one more unknown function than there are equations. The standard 
P/v approximation is simply to set the highest Legendre moment of ip(x,(i,t) equal to zero. 
In the case of the first four equations (corresponding to n = 0,1,2,3), one sets <f>4(x,t) = 
to get the following classic planar geometry time-dependent P3 equations: 



(x,t) + 



v dt 
v dt 



v dt 

Tx {t M 
d 



Ox 



dx V 5 

1 U( 



Q(x,t), 


(6a) 


0, 


(6b) 


0, 


(6c) 


0. 


(6d) 



From now on, to keep the discussion simple, we will work with the specific system of eqs. ([6]). 
However, everything which is done in the following can be generalized to more (or fewer) than 
four angular moments of eq. ([3]). 

Solving for the odd moments of eqs. Q, we divide eq. (6b) and eq. (6d) by S 4 (x) and 
obtain 



<f>i (a?, *) = --(/ + T)- 1 X{2(t)2{x, t) + Mx, t)), 



(j) 3 (x,t) 



(i+r)- 1 A'(30 2 (x,i)), 



(7a) 
(7b) 



where 



T 



1 







vT, t (x,i) dt 



and X 



1 







St(x, t) dx 



are two dimensionless operators. Introducing these expressions into the first and third of 
eqs. ([6]), we get 

llt( X,t ^~l~Hx~( I + 7 T 1 ' V ( 2< ^( X ' *) + Mx, t)) + E«(x, t)Mx, t) = Q{x, t), (9a) 



- — (x t)- — 
v dx ' dx 



~(/ + ry^Mx, t) + ~(i + T)- l x{2 ( t )2 (x, t) + Mx, t)) 

35 15 



+ E t (x,t)c/)2{x,t) = 0. (9b) 
To derive the SP3 equations, we formally replace the 1-D operator 

|(/ + n-'*-|(/ + r)-'^| do.) 
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by the 3-D operator: 



V • (J + T) 



-i 



1 



E t (x,t) 



-V 



(10b) 



9 ,-r ^ i 1 9 d /T n 1 9 9 ,, _ x i 1 9 . 

This step is purely ad hoc; yet it is exactly what is (or rather, was) done in the original (1960) 
derivation of the steady-state SPn equations [151 EH E! • Eqs. Q then become: 

-l^rfe *) + ^a{x, t)M%, t) = V ■ (J + T)- 1 --^— V(0o(x, t) + 20 2 (x, t)) + Q{x, t), 
v ot 61jt{x,t) 



Hal 



i_a^ 



•(x,t) + s t (x,t)02(x,t) = v- (/ + 7T 1 



i 



: E „,x,^(^ ofet) + y^ fet) ) ,llb) 



Remark 1. I£gs. |il) directly reduce to the standard steady-state SP3 equations when 

1 d 



T 



vY, t (x,t) dt 



0. 



Next, we rewrite eqs. (11) as a system of hyperbolic PDEs by defining 

1 



J (x,t) :=-(J + T) 



Z 2 fei) :=-(I + T) 



-1 



-1 



Then, eqs. (j 1 1 p can be written as: 

~{x,t) + V- J (x, t) + T ta (x,t)4> (x, t) = Q(x,t), 
'-(x,t) + V- J 2 (x,t) + Z t (x,t)(f>2(x, t) = 0, 



1 U<P0 , 

v dt 
ld<h. 
v dt 

' ' ^{x,t) + -V( </»ofet) + 2<A 2 fet) ) +St(x,t)J (x,t) = 0, 



3" 

1, 

3" 



v dt 

' ' ~ 2 "fe t) + ^Y (\<!>0{X, t) + y </> 2 (£. / ! ) +S,(£./)J 2 (>./.)= 0. 



(12a) 
(12b) 

(13a) 
(13b) 
(13c) 

(13d) 



These are time-dependent 3-D SP3 equations. This coupled system of first-order PDEs is 
easy to discretize, spatially and temporally. We note that in eqs. (13) all quantities (cross 
sections, source term and fluxes) can be functions of x and t. 

Remark 2. It is interesting to observe that the steady-state SPn equations are problematic 
in systems containing void regions (in which Sj(x,t) = 0), but eqs. (13) do not have an issue 
with voids. 
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Remark 3. Adding the partial time derivative to each of the steady-state SP3 equations in 
first-order form formally implies [31]: 



-— -(x,t) + V • <j) (x,t) +T, a (x > t)(j)o(x,t) = Q(x,t), 
v at — 1 

l^{x,t) + ^V- (20 1 (x,t) + 30 3 (x,t)) +£ t (x,t)0 2 (x,t) = 
1 dfa 1 

-^fet) + -Y(Mx,t) + 2fa(x,t)) + E t (x,t)^(x,t) = 0, 



v dt 



(x, t) + V ( -Mx, t) ) + E t (x, *) = 0, 



(14a) 
(14b) 
(14c) 
(14d) 



where the three-dimensional vectors i) and 4> 3 (x,t) are obtained by formally replacing 
the scalars (f>i(x,t) and <f>s{x,t) in eqs. (by by vectors. 



Although eqs. (14) differ from eqs. (TSy it is possible to transform them into each other 
by a similarity transformation acting on the PDE system of eqs. (IS). This transformation 
is defined by 



(15) 



Consequently, eqs. (13) and eqs. (14) are equivalent and their solutions are either identical or 
can easily be transformed into one another. 

The above remark confirms that numerical codes, based on time-dependent SP3 equations 
obtained by simply adding the time-derivative, in fact do approximate equations which are 
equivalent to those derived in this paper. However, since all existing codes are based on 
equations derived in a purely heuristic way, we now wish to provide a theoretical foundation 
for these equations, which are already being solved numerically. 



3 Formal Asymptotic Derivation 

To keep our discussion simple, we restrict the asymptotic analysis to the monoenergetic, 
3-D isotropically scattering particle transport in a homogeneous medium. However, a similar 
analysis might also be performed for anisotropic scattering, which is already presented in 
|26j for the steady-state equations. As in Section [2j we start our analysis with the linear 
Boltzmann equation: 

ldV 1 

-^(x,n,t) + n-v^(x,n,t) + z t v{x,n,t) = —(^ s Mx,t) + Q{x,t)), (i6a) 

V Ot 47T 

or, dividing by £t 

T* (x, n, t) + n ■ X}B(x, o, t) + * (x, 0, *) = t- ( c Mx, t) + Q ^ - ) , (i6b) 

4vr \ ^ J 
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then eq. (22) can be written as 



n,t) + MQ- x^(x, n, t) = — (f> (x, t), 

Air 



(24) 



or 



47T 



(25) 



Introducing eq. (25) into the definition of the current in eq. (18c): 



J(x,t) = / spz>(x,n,t)dn 

Js 2 

= — — f n(i + m n ■ x)~ Vofe t)dn, 

47T J S 2 



and then introducing this result into eq. (19), we end up with: 



1 

v dt 



-fe*) + T- I n- v( J + Mn-x)~ 1 Mx,t)dn + 'E a cf)o(x,t) 

47T J S 2 



(26) 
(27) 

Q(x,t). (28) 



This equation for <j>o(x,t) is formally exact, but the integral term is a complicated operator 
which needs further simplifications. We introduce a small, positive, dimensionless parameter 

e into eq. (28) such that the operator — V becomes small, i.e., 

St 



X 



(29) 



Remark 4. It should be emphasized that the only assumption for the following asymptotic 
analysis is that X^ = 0(e). Neither the time derivative is supposed to be small nor source 
terms are scaled. This is purely formal and is chosen to keep the framework as general as 
possible. To draw a line from the scaling considered in this paper to scalings from previous 
asymptotic SPjy derivations in literature, we list some of them which are acceptable for our 
asymptotics: 

• conventional scaling JT ffi [25 \ \27[ Pffi Here, the system is assumed to be scattering- 
dominated with its collision rate being much larger than its absorption rate. In this case, 
e is the ratio of the mean free path (which corresponds to T,^ 1 ) over a typical length scale 
for the solution. This scaling is the standard scaling, which has been used to gain the 
diffusion or steady-state SPn equations by performing an asymptotic analysis for e ~ 0. 



generalized conventional scaling [26]: Larsen introduces an alternate scaling which 
is physically consistent with the conventional scaling and additionally, includes free pa- 
rameters. Depending on the choice of these parameters either the standard or modified 
diffusion and SPjy equations are obtained. From the theoretical point of view, the latter 
equations are proven to increase the accuracy for deep penetration problems. 
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parabolic scaling JiOj/ : An asymptotic approach for the time- dependent SP3 equations 
is presented in which the time derivative, source term and absorption cross section 
are scaled by e 2 and the space derivative by e. This scaling can also be achieved by 
introducing 

(30a) 
(30b) 

(30c) 
(30d) 

where v, at, o~ a , q are of 0(1), into eq. and dividing by Ej. In addition to the physical 
assertions from the conventional scaling, it also requires that particles travel at high 
velocities. Combined with a high collision rate, low absorption rate, and small source 
terms, this scaling as a whole implies a slowly varying solution in space and an even 
smaller variation in time. 





V 


V 






£ 


s t 


_ a t 


e 




= £o- a , 


Q(x,t) 


= eq(x 



Next, we asymptotically expand the operator 

l = — / n • v(i + Mn- xy l dn. 

4vr J S 2 



(31) 



Due to the assertion in eq. (29) about the dimensionless spatial gradient X = e^V or any 
scaling which yields X_ = 0(e), we can expand the operator C in eq. (31 ) in a Neumann series. 

Thus, for e sufficiently small, we have that 

00 

C = X>irA>, (32) 

where 



n=0 



r 



— / 0- v(Mn-xy 

47T J S 2 



0{e r 



To achieve an 0(e 7 ) approximation we need the first seven C n : 
£ = 0, 

A = lY-il + T)- 1 * 
£2 = 0, 

£3 = ^ [Y • (/ + T)- l X] (I + 7T 1 [X-(I + T)- l X] 



C 4 



0, 
44 



£5 = ^ [Y • (/ + T)- 1 *] (/ + 7T 1 [X-(I + T^X] (I + 7T 1 [X-(I + T^X] 

U = 0, 

which are calculated in detail in the Appendix. 



(33) 

(34a) 
(34b) 
(34c) 
(34d) 
(34e) 
(34f) 
(34g) 
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Remark 5. It is important to emphasize one crucial assumption made in the derivation of 



eq. (34 d) and eq. (34J): Both are only exact for either homogeneous media or a system in 



which St depends only on one spatial variable. However, the rest of eqs. (34) are also exact 
for heterogeneous media. 



3.1 SPi Equations 

Ignoring terms of 0(e 3 ) in eq. (32) we obtain from eqs. (34) 

1 



V- {I + T)-'X + 0{e 6 ). 



Introducing this approximation for C in eq. (31) and eq. (28) we get 



- ^ (x, t) - • (I + TT^ofe t) + S a o (x, t) + 0(e 3 ) = Q(x, t). 
v at 3 



If we additionally define 



lo(x,t) ■■=-\{I + T)- 1 X(l )0 {x,t) 



and drop the error term, the above equations simplify to the SP\ equations: 
' )( ^ (x,t) + V • J (x,t) + T, a (x,t)(j)o(x,t) = Q(x,t), 



1 d J 



(35) 



v dt 

1 

v dt v "' ; 3 

Remark 6. Whereas the steady-state SP\ approximation is the standard diffusion equation 



(36a) 
(36b) 



which requires only one scalar-valued function (f>o, eqs. (36) are a system of two equations with 



the scalar variable 4>q and vector J . By simply adding the time- derivative to the steady-state 
SP± equation one might expect a parabolic time- dependent SP\ approximation which is also 
obtained in fld\j . However, above time- dependent SPi equations are hyperbolic. 



3.2 SP 2 Equations 



More accurate solutions can be gained systematically by taking higher order terms in eq. (32) 
into account. An asymptotically higher order approximation of 0(e 5 ) is given by: 



1<90 O , . 1 
v dt 3 



I + ^(I + T)- l L 



M*, *) + ^aMx, t) + 0(e 5 ) = Q(x, t) 



where, for the sake of notational simplicity, the operator L is defined as 

L := [X- (I + T)- 1 *] =0{e 2 ). 



We approximate the operator in square brackets in eq. (37) by 



and set 



I+-(I + T)- l L 



202 (x, t) + Mx,t) 



/--(/ + rr'L 



-1 



+ 0(e 4 ), 



0o (x, t), 



(37) 



(38) 



(39) 
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which can be rewritten to 
4 



15 



(/ + T)- 1 L(20 2 (x, t) + Mx, t)) = 2<f> 2 (x, t). 



(40) 



Hence, by applying the operator (I + T) on the left, we get 



- JU 2 fe t) + S^ 2 (x, t) = • (I + T)- l X{<t> Q {x, t) + 2</> 2 (x, t)). (41) 
v at 15 



Combining eq. (|37|) and eq. (41) as well as discarding the error term, we obtain the system 

'(x,t)+V - J (x,t) + Z a (x,t)<p (x,t) = Q(x,t), (42a) 
'■(x, t) + V • J 2 (x, t) + Et(x, t)0 2 fe t) = 0, (42b) 
. *) + (to(x, t) + 2(t> 2 (x, t)j + St(x, t) J (x, t) = 0, (42c) 

(42d) 



1 O^o 
1 a^y 2 

iaZo 

i^(x,t) + ^V( <>o(>./) +2<>- 2 (.r.l) ) +S,(£./)J 2 (.£./) =0. 



.X, 



3.3 SP3 Equations 



Having collected all operators, we approximate £ by truncating the series in eq. (32) at n = 6 



and introducing eqs. (34) into eq. (32): 

c = -\v-(i + ry 1 K 

4 



[V • (I + T)~ l X\ (I + 7T 1 [#•(/ + 7T 1 *] 
- [V-a + T)- 1 ^] ((/ + T)" 1 [^•(/ + T)- 1 £ Y]) 2 + 0(e 7 ) 

I v • (i + rr 1 * [1 + ^(1 + T)- 1 [X-(I + T)- l x] 



45 
44 

945 
1 

3 



We rewrite eq. (43d) to 



£ = - h t L {/ + ±(I + TY l L + i| ((J + T)- 1 ^ 2 } + 0(, 7 ) 

(/ + T)- 1 L|+0(e 7 ). 



4 

15 



(43a) 
(43b) 
(43c) 

(43d) 

(44a) 
(44b) 



Approximating the term in square brackets of eq. ( 44b ) like in eq. ( 38 ) , we conclude 
1. 



£ 



:£*£</ + 



15 



{I + T)- l L\+0{e 7 ) 



(45) 
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Using this approximation for C in eq. (31) and eq. (28), we get: 
1300, 



v dt 



(x, t) + T, a (f) (x, t) = Q(x, t) 



If we define 



(x, t) + 



2cf> 2 (x,t) :-- 



15 



(I + T)- 1 LM^t) \+0{e"' 



15 



(I + T)- 1 LMx,t), 



then 4>2 satisfies 



Mx,t) = —{I + T)- l L<j) Q {x,t). 
15 



Operating by (/ + T) on the last equation yields 
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fo{x,t) = —L(j) Q (x,t), 
15 



or 



IT +T)o 2 (.r.l) - L [ ^ Q {x,t) + ^M^t) 



which is equivalent to 
1 d 



1 



11 



v u dt 

Altogether, we get 
I30o 



L -(/>o{x,t) + —4>2(x,t) 



11 



+ S t 2 fe t) = -V • (I + 7T 1 — V -0 o fe t) + ^0 2 fe i) 



1 



(x, t) + S a o fe i) = V • (I + Ty^V (0 o (x, t) + 20 2 (x, t)) + Q(x, t), 
v dt Sht 



1 



1 



dt 



(x,t) + S t 2 (x,t) = V-^ + r)- 1 — V -0o(x,t) + — 2 (x,i) 



3E f 



11 



7 



(46) 

(47) 

(48) 

(49) 
(50) 

(51) 

(52a) 
(52b) 



In a homogeneous medium, eqs. (52) are identical to eqs. ( |ll[ ). This essentially proves the 
asymptotic derivation of the previously obtained eqs. (11). One can continue rewriting these 
equations to eqs. (13) in the same way as is done in Section [2] 



4 Numerical Results in 2D 

We perform 2-D simulations for diffusion, the Pn, and the SPn equations. The computations 
for the latter two approximations are done with a version of the code StaRMAp by Seibold and 
Frank |38j. The name StaRMAp stands for "Staggered grid Radiation Moment Approximation", 
which describes the key methodology of the approach. More specifically, it is a second order 
accurate finite difference method for linear hyperbolic balance laws of the form 

dtu + M x ■ d x u + My ■ d y u + C ■ u = q , (53) 
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where the matrices M x , M y , and G possess specific patterns of their nonzero entries, as 
described below. The numerical method is implemented in a concise Matlab code that the 
authors plan to make publicly available upon submission of the corresponding paper [38] . Let 
the components of the solution vector u be indexed by {1, 2, . . . , S}. The requirement on 
the nonzero entry patterns of M x , M y , and G is that the components of u can be distributed 
into four disjoint sets, according to {1, 2, . . . , S} = loo U Iio U loi U In, such that the following 
properties hold: 

(M x )i,j = 0V(i,j) i ((loo x /io) U (Iio x /oo) U (I i x In) U (In x I 01 )) , 
(M y ) id = V (i, j) £ ((loo x J i) U (Joi x loo) U (Iio x In) U (In x J 10 )) , (54) 
C id = V (i, j) £ ((loo x loo) U (Iio x Iio) U (loi x I i) U (In x In)) . 

With this distribution of the indices of the solution components, we consider the following 
four fully staggered sub-grids 

Goo = {(iAxjAy) \i,jeZ}, G w = {((i + \)Ax,jAy) \i,jeZ}, 
Goi = {(iAx, (j + I) Ay) G Z} , G n = {((i + I) Ax, (j + |)Ay) | i, j € Z} , 

and assign the components with indices in 1^ to the corresponding sub-grid Gki, where 
k,l G {0,1}. On these fully staggered grids, any spatial derivative is approximated by a 
simple central difference stencil: 

d x w(iAx,jAy) « i («;((« + l)Ax, jAy) - w((i - \)Ax,jAy)) 
d y w(iAx,jAy) « ^ (^(tAa;, (j + |)Ai/) - ^(iAx, (i - I) Ay)) e \7L . 

Hence, x-derivatives of components on G^e live on Gi_^ £, and y-derivatives of components 



on Gki live on Gk,i-e, where k,£ € {0, 1}. The nonzero entry patterns (54) guarantee that 
the distribution of the indices of u into the sets loo, Iio, Ioi 5 an d In is identical to the 
corresponding distribution of the indices of M x ■ d x u + M y ■ d y il + C ■ u. Both the classical 
Pjv equations as well as the here derived SPjy equations, possess precisely the nonzero entry 



patterns (54), that admit a solution on fully staggered grids. 



The time-derivative in (53) is resolved by bootstrapping. Having a time step At, we 
associate the components that live on Goo U Gn with the times To = {nAt \ n E Z}, and the 
components that live on Gio U Goi with the times T\ = {(n + \)At \ n G Z}. A full time step 
consists of two sub-steps: first, update information on the grid GioUGoi from time (re— \)At 
to (n + \)At, where information on Goo U Gn at the mid-time nAt is used; second, update 
information on the grid GooUGn from time nAt to (n+l)At, where information on GioUGoi 
at the mid-time (n + V)At is used. 

Specifically, the sub-step update rule is implemented as follows (here for Goo U Gn; the 
other sub-step works analogously). The terms r = q — M x ■ d x u — M y ■ d y u that come from 
Gio U Gqi at the mid-step time are considered constant over the sub-step. Thus, equation 



(53) becomes the ODE 

dtu + C-u = f (55) 
with r = const. In the special case that G is a diagonal matrix, G = diag(ci, . . . , eg), we can 



solve (55) from nAt to (n + l)Ai explicitly: 



u k (x, (n + l)At) = exp(-c k (x)At)u k (x, nAt) - i(l - exp(-c fc (x)At))r A . . (56) 
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Note that, if C is time-dependent, we evaluate it at the mid-step time, in which case (56) 



becomes a second-order accurate approximation to the true solution of (55). This second- 
order convergence was verified by the method of manufactured solutions [33] for the SPn 
equations of order N = 1,3. 

As Pjy as well as SPn solutions are spatially discontinuous at material interfaces for even 
iV [3D], we only present numerical calculations for odd N in the following. In some applica- 
tions the diffusion equation is used to calculate approximations to the Boltzmann equation. 
Concerning computational effort, this approach is comparably inexpensive for coarse to mod- 
erate spatial resolutions, in which case the resulting linear systems can be solved rapidly. 
Unfortunately, its accuracy is not satisfactory, as will be shown for example in Section |4.3[ 
For the sake of comparison, we solve the following time-dependent diffusion equation 



1 



0, V,t) = V • [D(x, y, t)V<t> (x, y, t)} - S a (z, y, t)(/> Q (x, y, t) + Q(x, y, t), (57) 
1 



D(x,y,t) 



(58) 



which can be interpreted as a low-order approximation to the Boltzmann equation |10| . 
We apply the second order central in space - Crank-Nicolson in time Finite Difference scheme 
to discretize above equation. 



Q(x, y,t) = ( 2 + sin(47rt)e 



-t/3 



4.1 Equivalence of SPn and Pn equations 

The detailed analysis of the (simplified) spherical harmonics approximation brought along cer- 
tain conditions under which the steady-state Pn and SPn equations are equivalent. McClar- 
ren describes some of them and demonstrates numerically this equivalence on a square of 
dimension L = 5 in a homogeneous medium with isotropic material coefficients and sources 
[31] . We slightly modify the inhomogeneous source from [31] and change it to the following 
time- dependent sinusoidal term 

1, (x,y) G [1.75,2.25] x [1.75,2.25], 
1, (%,y) 6 [2.75,3.25] x [1.5,2.5], 
1, (x,y) G [1.75,2.25] x [2.75,3.25], (59) 
1, (x,y) G [3.5,4.25] x [3.5,3.75], 
0, otherwise. 

All criteria for the equivalence were developed for the steady-state equations, and there is 
no obvious reason why they should also be true for the time-dependent case. We implement 
this test case with St = 1 and S a = 0.9 and investigate the temporal as well as steady- 
state behavior of various approximations. Periodic boundary conditions are enforced for all 
methods. 

This problem has an absorption coefficient which is nine times larger than the scattering 
coefficient. Due to the inhomogeneous source, the solution has additionally large spatial 
gradients. Consequently, an agreement between Pjy and SPn solutions cannot be justified by 
any of the asymptotic scalings mentioned in Section [3] 

Fig. [T] displays the scalar flux 4>q at t = 1 from several methods. Note that 4>$ is not in 
steady-state at this time. As the condition of small spatial derivatives, which underlie the 
diffusion approximation, is not met in this situation, the diffusion solution is inaccurate. The 
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(c) Diffusion (d) P39 



Figure 1: Approximations to (fro at tfi na i = 1, 250 x 250 discretization points. 

SPg and Pg solutions appear to be identical, which is also confirmed in Fig. [2] showing the 
scalar flux along x = 2 for (simplified) spherical harmonics of order N = 1,3,5,9. To the 
eye, there is again no difference between SPjv and P/v. Moreover, we observe a significant 
improvement from the SP\ to the SP3 approximation. Although SP5 still shows significant 
deviations, an order of N = 9 is sufficient for the simplified P/v approximation to be very 
close to the high-order P39 solution, which is considered as reference here. 

It turns out that all SPn and P/v solutions are indeed equivalent. To verify this behavior 
more precisely we calculate the difference in the L°°-norm 

max \4>Q PN {xi, yj) - (j)Q N (xi,yj)\ 

at several times for = 1,3,5,7,11 and observe that all differences are as small as the 
machine precision. Hence, not only do the P/v and SP/v equations have the same solutions 
but also the numerical discretizations of them, which is an advantageous property of the 
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0.5 1 1.5 2 2.5 3 3.5 4 4.5 

y 

Figure 2: 4>o at ^finai = 1 along x = 2, 250 x 250 discretization points: Pi (green dotted line), 
SPi (green crosses), P3 (purple dash-dot line), SP3 (purple triangles), P5 (blue solid line), 
SP5 (blue diamonds), P9 (red dashed line), SPg (red circles), P39 (black solid line) 

applied numerical scheme. Our numerical experiments indicate that the analysis, performed 
for the steady-state equations, to prove the equivalence of Pat and <SP/v equations in a general, 
homogeneous medium with isotropic cross sections and sources might be extended to the time- 
dependent equations including time-dependent sources. This is an issue of future work which 
goes beyond the purpose of this paper. 

Note that in this case of a 2-D problem the P/v equations have ( N + 1 )( N+2 ) unknowns, 
whereas the SPjv method consists only of 3( [^\ + 1) unknowns. Nevertheless, it is remarkable 
that the P/v and SP/v solutions in this example agree even for the time-dependent case. Hence, 
there exist situations where high-order SPn approximations yield very accurate solutions 
(here for N = 9), and the diffusion approximation gives unsatisfactory results. 

4.2 A Moving Rod 

The prediction of the behavior of nuclear reactors is essential for the design and safe operation 
of nuclear power plants. Apart from many other interactions, one challenge is to develop 
efficient and accurate techniques for the description of neutron distributions. 

We consider the linear transport Boltzmann equation ([!]) with slight modifications as 
before. Neglecting energy dependence and the coupling to precursors, approximative models 
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to the following equation on [—^,^] x [— ^, ^] are solved: 
1 <9* 

- -g^r 0> y, 0, + • W(z, y, n, t) + [E s + E/ + E 7 (x, y, t)]*(x, y, O, i) 



E s + vT,f 
4ir 



[ W(x,y,tf,t)dtf. (60) 
Js 2 



In particular, the initial condition for the scalar flux is set to 4>q(x, y, 0) = 10 6 and all 
other variables to zero. Periodic boundary conditions are enforced, and we use the following 
parameters: 

17 = 1, L = 2, i/ = 0.9, £/ = 2, E s = l. (61) 

In addition to the scattering process (taken into account by E s ), the above equation 
includes two relatively frequent interactions: T,f is the fission cross section describing the 
probability that a neutron will initiate a fission event when it collides with a nucleus; E 7 (x, y, t) 
is the capture cross section characterizing a capture event in which the nucleus gains a neutron. 
Consequently, the absorption cross section is the sum of both: 

E a (x, y, t) = E/ + E 7 (x, y, t). 

In a fission event, the target nucleus splits into two daughter nuclei and v is usually the mean 
number of fission neutrons that are released. However, to keep the test case simpler, we use 
< v < 1 as a parameter and set the capture cross section to 

E 7 (x,y,t) = (i/-l)S/ + ^ o (62) 

where the domain il^j is defined as 

= {( x , y) G M 2 : -0.5 < s < 0.2, -0.5 < y < 0, (63) 
-0.3 <x<0, 0<y< 0.6, (64) 
< x < 0.5, < y < 0.3}. (65) 

The reason for the complicated design of is that a simple rectangular shape would not 
produce spatial gradients which are large enough to distinguish between P3 and SP3 solutions. 

This definition of E 7 models an asymmetric rod with a cross section shown in Fig. [3^i, 
which is moved into or out of the moderator material in a way specified by the function s(t). 
We choose a sequence of three processes: 

• pushing the rod into the material over the time of AT, 

• keeping it in the moderator for the time of T max — 2AT and 

• pulling the rod out over the same time of AT, 
and, hence, define 



■t, < t < AT, 



AT "> 

s max , AT < t < T max - AT, (66) 



ma>; 

AT 



(T ma x j T max AT < t ^ Tj 
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X 



(a) Rod Geometry (b) P39 solution at t = 0.3: s max = 10, AT = 0.2 

Figure 3 

Tmax is the duration of the whole process and is set to T max = 0.6. The quantity s max 
determines how deep the rod is pushed into the moderator and can be associated with the 
maximum penetration depth. The ratio s max /AT models the speed of motion of the rod. For 
large values of s max /AT the rod is moved very quickly, and thus large time derivatives are 
generated. Due to the finite geometry of the rod, large gradients in space additionally occur. 

The main goal of this problem is to analyze the behavior of diffusion, SP/v and P/v solutions 
in 2D for large time and space gradients in a semi-realistic setting where the absorption cross 
section depends on all three variables x, y and t. 

Figure [4] displays a cut of the scalar flux distribution along y = for 251 x 251 discretization 
points. The P39 solution is considered as our reference, which can be seen in Fig. |3| When 
the control rod is moved into the moderator, neutrons are absorbed and hence, their number 
diminishes. If this is done at a small speed (compared to the velocity of the particles), 
neutrons will quickly flood the region close to the absorber. Therefore, the scalar flux (which 
is a quantity for the particle number) is smooth but still forms steep slopes in our setting 
(Fig. [4]) . Although the solution drastically decreases, both the SP3 and P3 approximations 
are close to the reference P39 result. Moreover, P3 is very close to SP3. On the other hand, 
the diffusion solution is inaccurate in regions where large spatial gradients are formed. For 
increasing time, after the rod has been pulled out of the system, the neutrons spread in the 
whole domain and their distribution flattens more and more until, in steady-state, (fro levels 
off to a constant value which is smaller than at the beginning of the process. Fig. [Ija shows 
the distribution at an advanced time, where P3 and SP3 are still accurate and the diffusion 
approximation is far off the reference solution. 

Our observations described above are consistent with the analysis in Section [3| The 
asymptotic analysis implies that the SP3 approximation is a higher-order correction to dif- 
fusion in cases where steep spatial slopes occur. As long as these slopes are not extremely 
large, SP3 and P3 approximations yield similar results. However, when a certain threshold is 
reached, the problem significantly leaves the asymptotic limit. As a consequence, differences 
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X 



Figure 5: Checkerboard test problem. Material coefficients: isotropic source (white), purely 
scattering E s = 1 = Sf (orange and white), purely absorbing £ a = 10 = S 4 (black). 

between SP3 and P3 become larger and SP3 might lose accuracy. 
4.3 Checkerboard 

We consider the checkerboard problem from [JJ: The domain is a square of size [0, 7] x [0, 7] 
where the majority of the region is purely scattering. In the middle of the lattice system, 
in the square [3,4] x [3,4], an isotropic source Q = 1 is continuously generating particles. 
Additionally, there are eleven small squares of side length one of purely absorbing spots in 
which S a = 10 = Tit- Fig. [5] illustrates the problem settings more precisely. Vacuum boundary 
conditions [1] are enforced and all initial quantities at t = are zero. We compare the scalar 
flux using different methods (including high-order Pn and SPn)- 

The diffusion solution in Fig. [6] gives a poor result. Particles are transported from the 
central source at a much higher speed to the boundaries of the domain. The solution is much 
smoother than it should be, and large slopes are washed out. 

Although the SP3 calculation shows large improvements upon the diffusion result, it is 
still too diffusive compared to the reference P39 (Fig. [6]) . In some purely absorbing regions 
away from the center, the P3 computation lacks particles because the approximation does not 
allow particle waves to travel at a high speed. Hence, depending on the desired purpose, P3 
is not evidently superior to SP3 in this case. Nevertheless, the particle beams between the 
corners of the absorbing regions are well resolved in all Pn solutions. 

For increasing N, the spherical harmonics solutions show better improvements than SPn- 
Although SP5 calculation yields visible changes, it still shows obvious differences from the 
reference. It is important to emphasize that SPn solutions of this problem cannot be expected 
to be highly accurate for several reasons: The underlying cross-sections are discontinuous and 
lead to a transport solution with steep slopes varying over at least seven orders of magnitude. 
Hence, the scaling parameter e from Section [3] cannot be small. Moreover, purely absorbing 
regions have a scattering ratio of c = ^ = which is significantly outside any asymptotic 
limit discussed in Section [3| And last, the total cross-section depends on both spatial 
variables x and y, which violates an assumption in our asymptotic derivation. 
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Figure 6: Checkerboard problem: Scalar flux <f)Q at t = 3.2 for from several approximations 
with 250 x 250 points. The values are plotted in log lo (0o) an d limited to seven orders of 
magnitude. 
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Figure 7: L -error between SPn and P39 or P/v and P39 for 1000 x 1000 spatial discretization 
points: e = 0.5 : P/v (black solid triangle line), e = 0.5 : SPn (black solid asterisk line), 
e = 0.1 : Pjv ( r ed solid cross line), e = 0.1 : SPn (red solid square line), e = 0.04 : P/v (blue 
solid circle line), s = 0.04 (blue solid plus line). 



However, realistic applications often include these kinds of geometries, in which cross- 
sections are (highly) varying and entail large spatial gradients in solutions. If SPn approxi- 
mations are computed for these problems, it will be of significant interest to know how much 
is gained by increasing order N. In other words, for growing N, how much does the er- 
ror between SPn and transport solution decrease? An answer to this question is necessary 
to determine the benefit of a high-order SPn calculation in comparison to the additional 
computational effort. 

To demonstrate the behavior of this error in a numerical way, SPn approximations are 
calculated for odd N, from N = 1 to N = 91. In order to make the numerical discretization 
error comparabely small we choose a very fine discretization grid of 1000 x 1000 points. 
Additionally, a scaling parameter e is introduced into all equations to achieve the scaling 
from eq. (1291). Approximative solutions are then computed for the scaled transport equation 



(19) where the velocity v is set to one. Again, P39 is assumed to be the reference solution. 

Fig. [7] shows the L 1 -error between SPn and P39 or P/y and P39 for different values of 
b = 0.04,0.1,0.5. The error is plotted against the system size (i.e., the number of equations 
to be solved) which is N + 1 for the SPn and (N + 1)(N + 2)/2 for the P/v approximation. 
This choice is made to give a reasonable comparison based on computational effort between 
Pjv and simplified Pn calculations. 

Indeed, both the Pn- and SP/v-error become smaller for decreasing values of e. More- 
over, for a fixed e and increasing order N, the error of the Pn and SPn method decreases 
monotonically. The P/v-error is even strictly decreasing to zero. In contrast, due to the dis- 
continuous property of the cross sections or a scattering ratio which is outside the asymptotic 
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limit, the S'P/v-error approaches a constant non-zero value and saturates at system sizes of 
roughly 10 - 20. 

For a fixed system size (i.e. the same number of equations to be solved), Fig. [7] shows 
iSP/v-errors which are smaller than the corresponding py-errors. Consequently, to achieve 
a desired L 1 -error above a minimum threshold the SPn approximation is computationally 
less costly than Pjy calculations. However, this is only valid for small system sizes and errors 
which are above the saturation limit of the SPn method. 

The asymptotic analysis in Section[3]is strictly valid for homogeneous media. Nevertheless, 
the SPn error function for a small e = 0.04 in Fig. [7] demonstrates that there is a significant 
error decrease for high order SPn solutions. Hence, this behavior shows numerically that 
SPn computations of large orders N > 5 can also be used in settings with discontinuous 
cross sections to achieve a specified error. 



5 Discussion 

We have presented an approach on how time-dependent SPn equations can be derived in 
an ad hoc way for general N, and have explicitely provided an asymptotic analysis for these 



equations up to N = 3. A crucial result is the system of SP3 equations from eqs. (13), a 



first-order system of hyperbolic PDEs. One feature of eqs. (13) is that in 1-D planar geome- 
try, they reduce exactly to the time-dependent 1-D planar-geometry P3 equations. Hence, in 
planar geometry, numerical solutions of the time-dependent SP3 and P3 equations are equiv- 
alent. Numerical calculations are therefore performed in two space dimensions in Section 
[4} Moreover, computations in Section 4.1 demonstrate that there are even cases outside the 



asymptotic limit where Pn and SPn are equivalent in 2D. 

Problems in heterogeneous media with cross sections of small scattering ratios c = |^ 
show the following behavior: The smaller the scaling parameter e becomes the closer are 
SPn computations to transport solutions. On the other hand, large e lead to large errors 



(Section 4.3). Solutions to problems which do not satisfy the homogeneity assumption, made 



in the asymptotic analysis, still show accurate solutions (Section 4.3). Furthermore, all com- 
parisons confirm that SP3 results are throughout superior to diffusion solutions, especially in 
cases where large gradients are formed. 



In the following, a few possible future tasks are discussed: 

1. The asymptotic analysis in this paper gives rise to time-dependent SP3 equations which 
differ from those developed in |10j . From a theoretical point of view, the main reason is 
the different scaling, which especially yields an additional e 2 in front of the time-derivative 
in [10]. Amongst others, one consequence of the different scaling is that the previous time- 
dependent SP3 equations do not reduce exactly to the time-dependent planar-geometry P3 
equations. Another difference is that the time-dependent SP3 theory in [10J contains an 
arbitrary constant a. Although all admissible values of a imply the same asymptotic order 
of accuracy, it could be still valuable to compare numerical solutions of those equations 
with numerical solutions of equations proposed in this publication. Central challenges are 
to explain different behaviors of these equations for certain problems and to verify whether 
computational results confirm the theoretical predictions. 

2. Recently, Larsen presented modified diffusion and SPn equations [26], which are designed 
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to be more accurate for deep penetration problems. The asymptotic analysis therein covers 
the steady-state, anisotropically scattering linear Boltzmann equation and is based on a 
different scaling. Can a similar procedure be applied to the time-dependent SPn equations 
considered in this paper? Answers to this question could provide a mathematical founda- 
tion for the time-dependent SPn equations in a much more general field of applications. 

3. It is difficult to calculate exact angular fluxes of the transport equation [J3] and obviously, 
even more difficult to compare them to approximations. Nevertheless, it is still possible 
to gain a different type of information about the analytic solution. One possibility to 
investigate the accuracy of angular moment approximations to the Boltzmann equation 
is to perform a moment analysis. The main idea of this method is to calculate angular 
and/or spatial moments of the angular flux and compare them to corresponding quantities 
in the moment equations. Consequently, one can draw conclusions about the accuracy of 
the method. Densmore and McClarren [8], e.g., compared the previous time-dependent 
simplified Pn methods from |10| to other approximations. A similar analysis would provide 
additional information about moments which are preserved by the time-dependent SPn 
solutions from this paper. 

4. McClarren suggests one important aspect about the SPn equations to be studied in the 
future |31j : What is the optimal order of SPn equations that determines a limit where 
your benefit in accuracy is still larger than the additional computational costs? Or in other 
words, up to which order N is it reasonable to use the SPn equations? No investigations 
concerning this issue have been performed up to now. However, the publicly available 
MATLAB code |38| used in this paper allows to perform Pn an d SPn calculations in 2D of 



very high order. In Section 4.3 the error behavior between SPn and the transport solution 
is studied for increasing order and different scaling parameters. The results indicate that 
up to a comparatively high order around N = 9 the SPn equations are superior to the Pn 
equations of same system size. These results could be a motivation for future analysis. 

5. As we only considered mono-energetic transport problems here, the extension to multi- 
group calculations along the lines of [29] will allow for the simulation of realistic applications 
and is left to future work. 



A Ad Hoc Derivation: Time-Dependent SPn equations of Ar- 
bitrary Order 

Although a mathematical foundation is only given for the SP3 equations in this paper, we 
also perform numerical computations with equations of higher order. Therefore, we briefly 
explain how time-dependent SPn equations of arbitrary order can be obtained. 

Starting from the 1-D time-dependent Pn equations in planar geometry, the following 
hyperbolic system of SPn equations can be derived by a similar procedure from Section [2] 

1 dcj) 

--^r(x,t) + V • J 2i (x,t) + Y, t (x,t)(j)2i(x,t) = 5 ifi (T, s 4> (x,t) + Q(x,t)), (67a) 
v ot 

1 dJ 

--ipfe*) + ¥_{ki4>2i+2(x,t) + k4>2i{x,t) + mi4>2i-2(x,t)} + T, t (x,t)J 2i (x,t), (67b) 
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where i = 0, 1, . 
by 



L^J and [ J is the floor function. The coefficients in eqs. (67) are given 



2i(2i - 1) 
(16z 2 — 1) ' 
32i 3 + 24i 2 



1 



(16i 2 - 1)(4» + 3)' 
4i 2 



mi 



[(16i 2 -l)' 
( 2(2i 2 + 3i + l) 
(4i + l)(4i + 3)' 



.0, 



N odd, 
iV even, 

2(t + 1) < JV, 
else. 



(68a) 
(68b) 

(68c) 



B Expanding Operators 

Here, we provide detailed algebraic manipulations for the simplification of the expansion 
operators 



c n = — [ n ■ v(M n ■ x) n dn, 

4vr J S 2 



(69) 



where 



Vljt Ul 2->t 4-7T 



]- I (-)dfi and M = (I + T)-\I -V). (70) 



If we denote the cosine of the polar angle by — 1 < fx < 1 and the azimuthal angle by 
< 4> < 2tt, then the unit vector 0_ = (fix, ^2^3) on the unit sphere is given by its scalar 
components 



(71) 
(72) 
(73) 



Q 2 = \A - c °s(</>), 
Q 3 = \J\ - /i 2 sin(</>). 

Much of the following analysis is based on results of the integral 



^71^22 • • • ^i n d£2 

s 2 

for all possible combinations of 1 < ii, . . . , i n < 3. Without a proof, we explicitly state some 
solutions which are needed later on: 



s 2 



J~2<£^ ^"^2 

is 2 



4lT 
47T 



15 



/ fi^Ojj . . . fii n dfi = for odd n. 

Js 2 



(74a) 
(74b) 
(74c) 



26 



We first focus on calculating C n for even n because above integral is zero for odd n: 

Cq = [ n-vdn = 0. 

4vr J S 2 

£2 = ^- I Q- Y(M n ■ x_f dn 

47T J S 2 

= — I n-v(Mn-x)(Mn-x)dn 

47T J S 2 

= t l / n-¥(Mn-x)(i + T)~ 1 (i-v)(n-x)dn. 

47T J S 2 



Since 



it follows 



C-2 



3 

7>(0 • x)^ fe t) = ^ V diMx, t) [ tii dn = 0, 

.- / 11- V(.MO •*)(/ + T) -1 ^ 

47T J52 
47T J52 

~ / o- vii + VHi^-^ii + V^in-x) 

4vr y 5 2 

- p (o • x) (j + 7T 1 (o • } dn 
- f + T)- 1 (n- x)(i + (n- x)dn 

:7r Js 2 



4ir j 52 
= 0, 

where last two equalities follow from the fact that 

[v (n-x) (i + TrHQ-x^M^t) 



(75) 

(76) 
(77) 
(78) 

(79) 

(80) 
(81) 

(82) 
(83) 
(84) 

(85) 



is independent of Q. Hence, only integrals of an odd number of are left which all vanish 



according to eq. ( 74c ) . Similarly, we obtain 



U = 0, 



(86) 
(87) 



and attend to odd-numbered operators C n : Note that eq. (79) yields 



= n-v(i + T)- 1 ^ r n-v. 
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Using eq. (74a), we obtain 

3 



4vr 7,52 4vr .^ =i J S 2 S t 

«1,«2 = 1 



J2 cZ^i 



3f", ,lv ' St 

^•(i + rr 1 ^. 



Next, making use of eq. (74c) we calculate 

£3 = ^/ n- v(A<n-^)(A<n-^)(A<n-^)dn 

47T J 5 2 

= f n-v(MQ- x)(mq- x)[(i + ry 1 (i -v)n - x]dn 

47T J 5 2 
47T J 5 2 

Am J S 2 

- — [ n- y(m n-x)(i + ry x v {n ■ x(i + r) _1 o • x] &i. 

47T J S 2 

According to eq. ( 74c ) , last integral simplifies to 

-— f n-v(Mn-x)(i + Ty 1 v{n-x(i + Ty 1 n-x}dn 

47T J S 2 

= -— f o- + ry 1 n- xm + ry 1 v \n- xu + ry 1 n- x\dn 

4vr y 5 2 

+ t- / ^■v{i + Ty l v{a-x{i + Ty 1 v{a-x{i + Ty 1 9 L -x}}dn 

4vr s v / 

=0 

= -— ! n-v[(i + Ty 1 n-x}{i + Ty 1 v\n-x{i + Ty 1 n-x}dn. 

4vr J S 2 

Similarly, we can eliminate the operator V in the first term of eq. (96) and obtain 

£3 = 7^ / + ry\i -v)^ - x][{i + ry 1 o L - x][{i + ry l n- x]dn 

4vr J S 2 

— — — f n-v[(i + ry 1 n- x](i + ry l v {n- x(i + ry l n- x}dn 

4tt J S 2 

= — f + ry 1 n- x][(i + ry 1 n- x][(i + ry 1 n- x]dn 

47T J S 2 

-— I + ry 1 n- x]{i + ry l v {n- x(i + ry l n- x\d£i 

47T J S 2 



S 2 
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We consider last two integral terms separately and expand the dot products therein: 



- f 0- + T)~ 1 n- x][(i + ry 1 n- x][(i + ry 1 n- xjdn 

^ .Is- 

l'l,... 1*4 = 1 ' 



(102) 



(/ + rr 1 — n i3 d i3 ) ( (/ + Ty l —n i4 d u ) dn (103) 



1 



-1 



1 



Si 



--itj^ (c+rr 1 ^,) (e+rr 1 ^) ((J+rr 1 ^) 



(104) 



E ^ (( J+r ) -1 ^^) (( /+r ^ 1 ^^) ( (I+rr1 ^) 
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(^1,12^*31*4 *^*1 i*3 ,*4 "I" ^*1 i*4^*2i*3 ) ' (10>5) 



where last equality follows from eq. (74b). Following the same arguments as in eqs. (89)-(90), 
the second integral term in eq. (101) can be rewritten to 



^ J • v [(/ + ry 1 9 L - x\{i + ty 1 v {n ■ x(i + ry 1 n-x} dn 

E fti (V+rr 1 ^) /+r ^ 1 ^) /+r ^ 1 ^ 1 



(106) 



. &ii ,12 $t 



*3i*4 - 



Altogether, we get 

3 



£3 = E ^((i + rr 1 ^) ^(j + r)- 1 ^) (V+^T 1 ^) 



45 



( 2(5^^2^43,14 "I - ,43 (5t a f i4 ~l~ 3<5j 1; j 4 (5 



*2,*3/ 



(107) 



Remark 7. Equation (101) is exact and holds for arbitrary geometries. In principle, eq. (101) 
can be used to derive equations which are asymptotically valid in heterogeneous media. How- 
ever, due to the complexity of eq. (101), one would probably lose the simple structure of the 
SPn equations as presented in this paper. To keep this main advantage of the SPn method 
further simplifications are needed and one possibility is shown in the following. 

Assuming that the medium is either 

• one dimensional (where the sum of 81 terms reduces to one single term) or 

• homogeneous (where Si is independent of the spatial variable x) , 
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we can rearrange the derivative operators in eq. (107) and conclude 



=^ [v • (i + rr 1 *] (/ + r)- 1 [# • (j + rr 1 *] . (ios) 

A similar analysis can also be performed for the next operator in the hierarchy. Since 
these expressions become even more involved and lengthy, we do not present them here. A 
road map to this analysis is given in [26.J where a detailed derivation is also presented for £3. 
Making the same assertions as above the operator simplifies to 

£5 = ^= [V • (J + T)- l X] (I + T)- 1 [X-(I + T)~ l X] (I + 7T 1 [X-{I + T)- 1 *} . (109) 
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